rice (Rice / Rician distribution)

The Rice (Rician) distribution is a continuous distribution on nonnegative real numbers. It arises as the distribution of the magnitude of a 2D Gaussian vector with a nonzero mean (equivalently, the amplitude of a complex Gaussian random variable with a deterministic “line-of-sight” component).

---

1) Title & Classification

Item	Value
Name	Rice (rice)
Type	Continuous
Support	\(r \in [0,\infty)\)
Parameters	noncentrality \(\nu\ge 0\), scale \(\sigma>0\)
Useful reparam.	\(K = \nu^2/(2\sigma^2)\) (Rician \(K\)-factor)

What you’ll be able to do after this notebook

• connect Rice to a shifted 2D Gaussian and to Rayleigh / noncentral \(\chi\) distributions

• write the PDF and CDF in standard special-function notation

• compute mean/variance/skew/kurtosis from raw moments

• derive the PDF (polar change of variables) and the log-likelihood

• sample from Rice using NumPy only

• use scipy.stats.rice for evaluation, simulation, and fit

import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio

from scipy import stats, special
from scipy.optimize import minimize

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=6, suppress=True)

rng = np.random.default_rng(42)

2) Intuition & Motivation

What this distribution models

A canonical construction:

• Let \((X,Y)\) be independent with $\(X \sim \mathcal{N}(\nu,\sigma^2),\qquad Y \sim \mathcal{N}(0,\sigma^2).\)$

• Define the amplitude (radius) $\(R = \sqrt{X^2 + Y^2}.\)$

Then \(R\) has a Rice distribution with parameters \((\nu,\sigma)\).

Intuition:

• When \(\nu=0\), the mean is at the origin, so the radius is the familiar Rayleigh distribution.

• When \(\nu>0\), the Gaussian cloud is shifted away from the origin; radii concentrate around \(\approx \nu\).

• The dimensionless ratio \(K = \nu^2/(2\sigma^2)\) (often called the Rician \(K\)-factor) measures how strong the deterministic component is relative to the random scatter.

Typical real-world use cases

• Wireless communications: Rician fading amplitudes (line-of-sight + multipath scatter).

• Radar / sonar: return amplitude with a deterministic component plus noise.

• MRI magnitude images: magnitude of complex Gaussian noise leads to Rician-like likelihoods.

• Signal processing: envelope of a sinusoid in additive Gaussian noise.

Relations to other distributions

• Rayleigh: \(\nu=0\) gives Rayleigh\((\sigma)\).

• Noncentral \(\chi\): Rice is the noncentral-\(\chi\) distribution with degrees of freedom \(k=2\).

• Noncentral \(\chi^2\): \(R^2/\sigma^2 \sim \chi'^2(2,\lambda)\) with noncentrality \(\lambda=(\nu/\sigma)^2\).

• Approx. Normal for large \(K\): when \(\nu\gg\sigma\), \(R\) is close to \(\mathcal{N}(\nu,\sigma^2)\) truncated at \(0\) (the truncation becomes negligible).

3) Formal Definition

We use the parameterization \((\nu,\sigma)\) with \(\nu\ge 0\) and \(\sigma>0\).

PDF

For \(r\ge 0\):

\[ f(r\mid \nu,\sigma) = \frac{r}{\sigma^2}\exp\!\left(-\frac{r^2+\nu^2}{2\sigma^2}\right) I_0\!\left(\frac{r\nu}{\sigma^2}\right), \]

and \(f(r)=0\) for \(r<0\). Here \(I_0(\cdot)\) is the modified Bessel function of the first kind (order 0).

A useful integral identity (often used in derivations) is:

\[ I_0(t) = \frac{1}{\pi}\int_0^{\pi} e^{t\cos\theta}\,d\theta. \]

CDF

A standard special-function expression uses the Marcum \(Q\)-function of order 1:

\[ F(r\mid \nu,\sigma) = 1 - Q_1\!\left(\frac{\nu}{\sigma},\frac{r}{\sigma}\right),\qquad r\ge 0. \]

One definition of \(Q_1(a,b)\) is

\[ Q_1(a,b) = \int_b^{\infty} x\,\exp\!\left(-\frac{x^2+a^2}{2}\right) I_0(ax)\,dx. \]

SciPy provides numerically stable evaluation of the CDF via scipy.stats.rice.cdf.

SciPy parameter mapping

SciPy’s Rice distribution is scipy.stats.rice with shape parameter \(b=\nu/\sigma\):

\[ R \sim \text{Rice}(\nu,\sigma) \quad\Longleftrightarrow\quad \texttt{stats.rice(b=nu/sigma, loc=0, scale=sigma)}. \]

def rice_logpdf(r, nu, sigma):
    '''Log-PDF of Rice(nu, sigma) for r in [0, inf).

    Uses exponentially-scaled Bessel i0e for numerical stability.
    '''

    r = np.asarray(r, dtype=float)
    nu = float(nu)
    sigma = float(sigma)

    if nu < 0 or sigma <= 0:
        raise ValueError("nu must be >= 0 and sigma must be > 0")

    logpdf = np.full_like(r, -np.inf, dtype=float)
    mask = r > 0
    if not np.any(mask):
        return logpdf

    z = (r[mask] * nu) / (sigma**2)

    # I0(z) grows ~ exp(z)/sqrt(2*pi*z), so work with i0e(z) = exp(-z) I0(z)
    log_i0 = np.log(special.i0e(z)) + z

    logpdf[mask] = (
        np.log(r[mask])
        - 2.0 * np.log(sigma)
        - (r[mask] ** 2 + nu**2) / (2.0 * sigma**2)
        + log_i0
    )
    return logpdf


def rice_pdf(r, nu, sigma):
    return np.exp(rice_logpdf(r, nu, sigma))


def rice_rvs_numpy(nu, sigma, size=1, rng=None):
    '''Sample from Rice(nu, sigma) using NumPy only.

    Construction: R = sqrt(X^2 + Y^2) where X~N(nu, sigma^2), Y~N(0, sigma^2).
    '''

    if rng is None:
        rng = np.random.default_rng()

    nu = float(nu)
    sigma = float(sigma)
    if nu < 0 or sigma <= 0:
        raise ValueError("nu must be >= 0 and sigma must be > 0")

    x = rng.normal(loc=nu, scale=sigma, size=size)
    y = rng.normal(loc=0.0, scale=sigma, size=size)
    return np.sqrt(x * x + y * y)


def ecdf(samples):
    x = np.sort(np.asarray(samples, dtype=float))
    y = np.arange(1, x.size + 1) / x.size
    return x, y

4) Moments & Properties

A convenient organizing principle is raw moments

\[m_n = \mathbb{E}[R^n].\]

Let \(a = \nu^2/(2\sigma^2)\).

Raw moments

For \(n>-2\) (in particular for \(n=1,2,3,4\)):

\[ m_n = \sigma^n\,2^{n/2}\,\Gamma\!\left(1+\frac{n}{2}\right) \,{}_1F_1\!\left(-\frac{n}{2};\,1;\,-a\right), \]

where \({}_1F_1\) is the confluent hypergeometric function.

From these:

• Mean: \(\mu = m_1\)

• Variance: \(\mathrm{Var}(R)=m_2-m_1^2\)

A commonly used explicit expression for the mean (in Bessel functions) is

\[ \mathbb{E}[R] = \sigma\sqrt{\frac{\pi}{2}}\,e^{-a/2} \left[(1+a)I_0\!\left(\frac{a}{2}\right)+a\,I_1\!\left(\frac{a}{2}\right)\right]. \]

A very simple second raw moment comes from the 2D Gaussian construction:

\[ \mathbb{E}[R^2] = 2\sigma^2 + \nu^2. \]

Therefore:

\[ \mathrm{Var}(R) = 2\sigma^2 + \nu^2 - \big(\mathbb{E}[R]\big)^2. \]

Skewness and kurtosis

From raw moments \(m_1,\dots,m_4\):

[ \mu_2 = m_2 - m_1^2,\quad \mu_3 = m_3 - 3 m_1 m_2 + 2 m_1^3,\quad \mu_4 = m_4 - 4 m_1 m_3 + 6 m_1^2 m_2 - 3 m_1^4. ]

[ \text{skewness} = \frac{\mu_3}{\mu_2^{3/2}},\qquad \text{kurtosis} = \frac{\mu_4}{\mu_2^{2}},\qquad \text{excess kurtosis} = \frac{\mu_4}{\mu_2^{2}} - 3. ]

SciPy can compute \((\mu,\mathrm{Var},\text{skew},\text{kurtosis})\) via stats.rice.stats(..., moments="mvsk").

MGF / characteristic function

The MGF of \(R\) exists for all real \(t\) (the tail is essentially Gaussian), but it does not simplify to elementary functions:

\[ M_R(t)=\mathbb{E}[e^{tR}] = \int_0^{\infty} e^{tr} f(r)\,dr. \]

A useful closed form exists for the squared amplitude \(R^2\):

• \(R^2 = X^2+Y^2\) with \(X\sim\mathcal{N}(\nu,\sigma^2)\), \(Y\sim\mathcal{N}(0,\sigma^2)\).

• Equivalently, \(R^2/\sigma^2 \sim \chi'^2(2,\lambda)\) with \(\lambda=(\nu/\sigma)^2\).

For \(t<1/(2\sigma^2)\):

\[ M_{R^2}(t) = \frac{1}{1-2\sigma^2 t}\,\exp\!\left(\frac{\nu^2 t}{1-2\sigma^2 t}\right). \]

The characteristic function follows by \(t\mapsto it\).

Entropy

The differential entropy is

\[ H(R) = -\int_0^{\infty} f(r)\log f(r)\,dr = -\mathbb{E}[\log f(R)]. \]

There is no simple closed form in general; it is commonly evaluated numerically (quadrature or Monte Carlo).

def rice_raw_moment(n, nu, sigma):
    '''Raw moment m_n = E[R^n] for Rice(nu, sigma), using 1F1 representation.'''
    nu = float(nu)
    sigma = float(sigma)
    if nu < 0 or sigma <= 0:
        raise ValueError("nu must be >= 0 and sigma must be > 0")

    a = nu**2 / (2.0 * sigma**2)
    return (sigma**n) * (2.0 ** (n / 2.0)) * special.gamma(1.0 + n / 2.0) * special.hyp1f1(
        -n / 2.0, 1.0, -a
    )


def rice_mean(nu, sigma):
    '''Mean of Rice(nu, sigma), computed in a numerically stable Bessel form.'''
    nu = float(nu)
    sigma = float(sigma)
    if nu < 0 or sigma <= 0:
        raise ValueError("nu must be >= 0 and sigma must be > 0")

    a = nu**2 / (2.0 * sigma**2)
    z = a / 2.0

    # mean = sigma*sqrt(pi/2)*exp(-z)*[(1+a)I0(z) + a I1(z)]
    # use i0e(z)=exp(-z)I0(z), i1e(z)=exp(-z)I1(z)
    return sigma * np.sqrt(np.pi / 2.0) * ((1.0 + a) * special.i0e(z) + a * special.i1e(z))


def rice_central_moments(nu, sigma):
    m1 = rice_raw_moment(1, nu, sigma)
    m2 = rice_raw_moment(2, nu, sigma)
    m3 = rice_raw_moment(3, nu, sigma)
    m4 = rice_raw_moment(4, nu, sigma)

    mu2 = m2 - m1**2
    mu3 = m3 - 3 * m1 * m2 + 2 * m1**3
    mu4 = m4 - 4 * m1 * m3 + 6 * m1**2 * m2 - 3 * m1**4
    return m1, mu2, mu3, mu4


def rice_skew_kurtosis(nu, sigma):
    mean, mu2, mu3, mu4 = rice_central_moments(nu, sigma)
    skew = mu3 / (mu2 ** 1.5)
    kurt = mu4 / (mu2**2)
    return mean, mu2, skew, kurt


# Quick numeric check vs SciPy
nu0, sigma0 = 2.0, 1.3
b0 = nu0 / sigma0

mean_ours, var_ours, skew_ours, kurt_ours = rice_skew_kurtosis(nu0, sigma0)
mean_scipy, var_scipy, skew_scipy, kurt_scipy = stats.rice(b=b0, scale=sigma0).stats(
    moments="mvsk"
)

(mean_ours, mean_scipy), (var_ours, var_scipy), (skew_ours, skew_scipy), (kurt_ours, kurt_scipy)

((2.474457072384077, 2.4744570723840766),
 (1.2570621969284224, 1.2570621969284217),
 (0.34234947509082275, 0.3423494750908338),
 (2.835812529966139, -0.1641874700338617))

5) Parameter Interpretation

Meaning of the parameters

• \(\nu\) (noncentrality / deterministic component)

  • Geometrically: distance of the 2D Gaussian mean from the origin.

  • In communications: proportional to the strength of the line-of-sight component.

• \(\sigma\) (scale / scatter per quadrature)

  • Standard deviation of each Gaussian component \(X\) and \(Y\).

  • Controls how “spread out” the cloud is and therefore how variable the amplitude is.

A common dimensionless summary is the \(K\)-factor

\[ K = \frac{\nu^2}{2\sigma^2}, \]

interpretable as a (signal power)/(scatter power) ratio.

Shape changes

• Fix \(\sigma\) and increase \(\nu\):

  • density shifts to the right, mass near \(0\) decreases,

  • skewness decreases, and the distribution becomes more “Gaussian-like”.

• Fix \(\nu\) and increase \(\sigma\):

  • density spreads out,

  • the mode moves left and the distribution becomes more right-skewed.

# How the PDF changes with nu (sigma fixed)

sigma_viz = 1.0
nus = [0.0, 0.5, 1.5, 3.0]

x = np.linspace(0, 8, 600)
fig = go.Figure()
for nu in nus:
    fig.add_trace(
        go.Scatter(x=x, y=rice_pdf(x, nu=nu, sigma=sigma_viz), mode="lines", name=f"ν={nu:g}")
    )

fig.update_layout(
    title=f"Rice PDF for varying ν (σ={sigma_viz:g} fixed)",
    xaxis_title="r",
    yaxis_title="density",
)
fig.show()

6) Derivations

PDF (polar change of variables)

Start from a shifted isotropic Gaussian:

\[\begin{split} (X,Y) \sim \mathcal{N}\left(\begin{bmatrix}\nu\\0\end{bmatrix},\;\sigma^2 I_2\right). \end{split}\]

The joint density is

\[ f_{X,Y}(x,y)=\frac{1}{2\pi\sigma^2}\exp\!\left(-\frac{(x-\nu)^2+y^2}{2\sigma^2}\right). \]

Transform to polar coordinates \(x=r\cos\theta\), \(y=r\sin\theta\) with Jacobian \(|J|=r\):

\[ f_{R,\Theta}(r,\theta) = \frac{r}{2\pi\sigma^2}\exp\!\left(-\frac{r^2+\nu^2-2r\nu\cos\theta}{2\sigma^2}\right). \]

Now integrate out the angle:

(26)\[\begin{align} f_R(r) &= \int_0^{2\pi} f_{R,\Theta}(r,\theta)\,d\theta\\ &= \frac{r}{\sigma^2}\exp\!\left(-\frac{r^2+\nu^2}{2\sigma^2}\right) \underbrace{\frac{1}{2\pi}\int_0^{2\pi}\exp\!\left(\frac{r\nu}{\sigma^2}\cos\theta\right)\,d\theta}_{=\;I_0\left(r\nu/\sigma^2\right)}. \end{align}\]

This yields the Rice PDF.

Expectation and variance

The squared radius has a simple form:

\[ R^2 = X^2+Y^2. \]

Using \(X=\nu+\sigma Z_1\), \(Y=\sigma Z_2\) with \(Z_1,Z_2\sim\mathcal{N}(0,1)\) independent:

(27)\[\begin{align} \mathbb{E}[R^2] &= \mathbb{E}[(\nu+\sigma Z_1)^2] + \mathbb{E}[(\sigma Z_2)^2]\\ &= \nu^2 + \sigma^2\mathbb{E}[Z_1^2] + \sigma^2\mathbb{E}[Z_2^2]\\ &= \nu^2 + 2\sigma^2. \end{align}\]

Then

\[ \mathrm{Var}(R) = \mathbb{E}[R^2] - (\mathbb{E}[R])^2 = 2\sigma^2 + \nu^2 - (\mathbb{E}[R])^2. \]

Likelihood (i.i.d. sample)

For data \(r_1,\dots,r_n\) (all \(\ge 0\)), the log-likelihood is

\[ \ell(\nu,\sigma) = \sum_{i=1}^n\Bigg[ \log r_i - 2\log\sigma -\frac{r_i^2+\nu^2}{2\sigma^2} +\log I_0\!\left(\frac{r_i\nu}{\sigma^2}\right) \Bigg]. \]

There is no closed-form MLE for \((\nu,\sigma)\) in general; numerical optimization is standard.

def rice_negloglik(params, data):
    nu, sigma = float(params[0]), float(params[1])
    if nu < 0 or sigma <= 0:
        return np.inf
    return -np.sum(rice_logpdf(data, nu, sigma))


# Synthetic data + MLE (numerical)
nu_true, sigma_true = 2.0, 1.0
n = 300

data = rice_rvs_numpy(nu_true, sigma_true, size=n, rng=rng)

# Start from SciPy fit as a reasonable initial guess (loc fixed at 0)
b_hat, loc_hat, scale_hat = stats.rice.fit(data, floc=0)
nu_init, sigma_init = b_hat * scale_hat, scale_hat

res = minimize(
    rice_negloglik,
    x0=np.array([nu_init, sigma_init]),
    args=(data,),
    method="L-BFGS-B",
    bounds=[(0.0, None), (1e-12, None)],
)

nu_mle, sigma_mle = res.x
(nu_true, sigma_true), (nu_mle, sigma_mle), res.success

((2.0, 1.0), (1.9683706333582578, 0.9633278327527912), True)

7) Sampling & Simulation

NumPy-only algorithm

Use the defining construction:

1. Draw \(X\sim\mathcal{N}(\nu,\sigma^2)\).

2. Draw \(Y\sim\mathcal{N}(0,\sigma^2)\) independent.

3. Return \(R = \sqrt{X^2+Y^2}\).

Why this works: \(R\) is the radius of a 2D Gaussian centered at \((\nu,0)\).

This approach is fast, vectorized, and avoids special-function inversion.

# Smoke test: Monte Carlo mean/variance vs theory
nu_test, sigma_test = 1.5, 0.8
s = rice_rvs_numpy(nu_test, sigma_test, size=200_000, rng=rng)

mean_mc = s.mean()
var_mc = s.var()
mean_th = rice_mean(nu_test, sigma_test)
var_th = rice_raw_moment(2, nu_test, sigma_test) - mean_th**2

(mean_mc, mean_th), (var_mc, var_th)

((1.7337792125752307, 1.7345364036221205),
 (0.5240102238815233, 0.5213834645096407))

8) Visualization

We’ll compare:

• PDF vs Monte Carlo histogram

• CDF vs empirical CDF

We’ll compute the theoretical CDF using scipy.stats.rice (stable implementation).

nu_viz, sigma_viz = 2.0, 1.0
b_viz = nu_viz / sigma_viz

dist = stats.rice(b=b_viz, loc=0, scale=sigma_viz)

samples = rice_rvs_numpy(nu_viz, sigma_viz, size=80_000, rng=rng)

x_max = float(np.quantile(samples, 0.995))
x_grid = np.linspace(0, x_max, 600)

# PDF + histogram
fig = go.Figure()
fig.add_trace(
    go.Histogram(
        x=samples,
        nbinsx=80,
        histnorm="probability density",
        name="Monte Carlo samples",
        opacity=0.55,
    )
)
fig.add_trace(
    go.Scatter(
        x=x_grid,
        y=rice_pdf(x_grid, nu=nu_viz, sigma=sigma_viz),
        mode="lines",
        name="Theoretical PDF (NumPy+SciPy special)",
        line=dict(width=3),
    )
)
fig.update_layout(
    title=f"Rice(ν={nu_viz:g}, σ={sigma_viz:g}): histogram vs PDF",
    xaxis_title="r",
    yaxis_title="density",
    bargap=0.02,
)
fig.show()

# CDF + empirical CDF
xs, ys = ecdf(samples)

fig = go.Figure()
fig.add_trace(
    go.Scatter(
        x=x_grid,
        y=dist.cdf(x_grid),
        mode="lines",
        name="Theoretical CDF (SciPy)",
        line=dict(width=3),
    )
)
fig.add_trace(
    go.Scatter(
        x=xs[::120],
        y=ys[::120],
        mode="markers",
        name="Empirical CDF (subsampled)",
        marker=dict(size=5),
    )
)
fig.update_layout(
    title=f"Rice(ν={nu_viz:g}, σ={sigma_viz:g}): empirical CDF vs CDF",
    xaxis_title="r",
    yaxis_title="CDF",
)
fig.show()

9) SciPy Integration (scipy.stats.rice)

Key points:

• Shape parameter: b = nu / sigma (dimensionless).

• Scale parameter: scale = sigma.

• Location shift: loc (use loc=0 for the standard Rice support \([0,\infty)\)).

Common methods:

• pdf, logpdf, cdf, ppf

• rvs (random sampling)

• fit (MLE; consider fixing loc=0)

nu_true, sigma_true = 2.0, 1.0
b_true = nu_true / sigma_true

dist = stats.rice(b=b_true, loc=0, scale=sigma_true)

x = np.linspace(0, 8, 6)

pdf = dist.pdf(x)
cdf = dist.cdf(x)

samples_scipy = dist.rvs(size=5000, random_state=rng)

# Fit (fix loc=0 to match the usual Rice support)
b_hat, loc_hat, scale_hat = stats.rice.fit(samples_scipy, floc=0)
nu_hat = b_hat * scale_hat

(pdf, cdf), (nu_hat, scale_hat)

((array([0.      , 0.346012, 0.250942, 0.012433, 0.000045, 0.      ]),
  array([0.      , 0.242633, 0.841085, 0.995869, 0.99999 , 1.      ])),
 (1.9749664963783742, 0.9964787814644809))

10) Statistical Use Cases

A) Hypothesis testing

A common test in communications / detection:

• \(H_0\): Rayleigh fading (\(\nu=0\))

• \(H_1\): Rician fading (\(\nu>0\))

A likelihood ratio test compares the best fit under \(\nu=0\) (Rayleigh) versus the best fit under \(\nu\ge 0\) (Rice).

B) Bayesian modeling

Rice likelihoods appear when you observe a magnitude but the latent signal is complex Gaussian. Typical Bayesian ingredients:

• priors on \((\nu,\sigma)\) (or on \(K\) and \(\sigma\))

• hierarchical models: different groups/locations share hyperparameters

• inference via MCMC/VI since the posterior is non-conjugate

C) Generative modeling

Rice is a natural emission distribution for nonnegative amplitudes generated from an underlying complex Gaussian latent:

• sample latent complex value (mean + Gaussian noise)

• emit its magnitude

This shows up in fading simulators and in synthetic MRI magnitude data generation.

# Simple likelihood-ratio statistic: Rayleigh (nu=0) vs Rice (nu>=0)

def rayleigh_sigma_mle(data):
    data = np.asarray(data, dtype=float)
    return np.sqrt(np.mean(data**2) / 2.0)


def loglik_rice(data, nu, sigma):
    return float(np.sum(rice_logpdf(data, nu, sigma)))


def loglik_rayleigh(data, sigma):
    return loglik_rice(data, nu=0.0, sigma=sigma)


# Simulate a dataset under H1
nu_alt, sigma_alt = 1.6, 1.0
x = rice_rvs_numpy(nu_alt, sigma_alt, size=250, rng=rng)

# Fit under H0 (Rayleigh) and H1 (Rice)
sigma0 = rayleigh_sigma_mle(x)
ll0 = loglik_rayleigh(x, sigma0)

b1, _, s1 = stats.rice.fit(x, floc=0)
nu1 = b1 * s1
ll1 = loglik_rice(x, nu1, s1)

llr = 2.0 * (ll1 - ll0)
(sigma0, (nu1, s1), llr)

(1.4926947250781553,
 (1.5823631882104543, 0.9880517248571086),
 12.58952187373336)

11) Pitfalls

• Invalid parameters:

  • Require \(\nu\ge 0\) and \(\sigma>0\) (SciPy: b>=0, scale>0).

• Parameterization mismatch:

  • Many references use \((\nu,\sigma)\); SciPy uses b=nu/sigma plus scale=sigma.

  • Some domains use \(K=\nu^2/(2\sigma^2)\); make sure you track units.

• loc in SciPy:

  • stats.rice.fit may introduce a nonzero loc unless constrained.

  • For the canonical Rice model, it’s common to fix loc=0.

• Numerical overflow in \(I_0\):

  • The Bessel function \(I_0(z)\) grows like \(\exp(z)\).

  • Prefer log-space computations and the scaled function scipy.special.i0e.

• Flat/ill-conditioned likelihoods:

  • When \(\nu/\sigma\) is very small, Rice is close to Rayleigh.

  • When \(\nu/\sigma\) is very large, \(R\) is close to normal around \(\nu\).

  • Both regimes can make numerical fitting sensitive to initialization.

# Numerical stability demo: naive PDF can overflow for large nu/sigma

nu_big, sigma_small = 30.0, 1.0
r_grid = np.linspace(0, 50, 6)

# Stable logpdf
lp = rice_logpdf(r_grid, nu_big, sigma_small)

# Naive pdf (may overflow internally via I0)
# Wrap in errstate: the point is to show unstable behavior without noisy warnings.
# NOTE: this is for demonstration; prefer rice_pdf / rice_logpdf above.
with np.errstate(over="ignore", under="ignore", invalid="ignore"):
    naive = (
    (r_grid / sigma_small**2)
    * np.exp(-(r_grid**2 + nu_big**2) / (2 * sigma_small**2))
    * special.i0((r_grid * nu_big) / (sigma_small**2))
)

lp, naive

(array([       -inf, -201.467827,  -51.121463,   -0.9188  ,  -50.774993,
        -200.663442]),
 array([ 0.,  0.,  0., nan, nan, nan]))

12) Summary

• Rice is a continuous distribution on \([0,\infty)\) modeling the magnitude of a shifted 2D Gaussian.

• Parameters: \(\nu\) (deterministic/LOS component) and \(\sigma\) (scatter scale); \(K=\nu^2/(2\sigma^2)\) is a common dimensionless summary.

• PDF involves a modified Bessel function \(I_0\); CDF can be written using the Marcum \(Q\)-function.

• Raw moments have a compact form via \({}_1F_1\); \(\mathbb{E}[R^2]=2\sigma^2+\nu^2\) gives a simple variance identity.

• Sampling is easy with NumPy: draw \((X,Y)\) as Gaussians and take \(\sqrt{X^2+Y^2}\).

References

• SciPy documentation: scipy.stats.rice

• M. K. Simon and M.-S. Alouini, Digital Communication over Fading Channels (Rician fading)

• Classic special-function relations for Rice / noncentral-\(\chi\) distributions